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We consider the singular value statistics of products of independent 
random matrices. In particular we compute the corresponding averages 
of products of characteristic polynomials. To this aim we apply the pro¬ 
jection formula recently introduced for chiral random matrix ensembles 
which serves as a short cut of the supersymmetry method. The projection 
formula enables us to study the local statistics where free probability cur¬ 
rently fails. To illustrate the projection formula and underline the power 
of our approach we calculate the hard edge scaling limit of the Meijer G- 
ensembles comprising the Wishart-Laguerre (chiral Gaussian), the Jacobi 
(truncated orthogonal, unitary or unitary symplectic) and the Cauchy- 
Lorentz (heavy tail) random matrix ensembles. All calculations are done 
for real, complex, and quaternion matrices in a unifying way. In the case 
of real and quaternion matrices the results are completely new and point 
to the universality of the hard edge scaling limit for a product of these ma¬ 
trices, too. Moreover we identify the non-linear er-models corresponding 
to product matrices. 

PACS numbers: 02.10.Yn, 02.50.Sk, 05.40.-a 

1. Introduction 

Sums and products of random matrices are the simplest generalization 
of random matrix theory (RMT) to introduce some kind of dimension. 
Sums of random matrices can be understood as a convolution and regu¬ 
larly appear in the field of Dyson’s Brownian motion pQ. Product matrices 
are versatile as well. Applications of them can be found in mesoscopic 
physics mm, qcd m, and wireless telecommunication urn- In the past 
years a lot of progress was made on products of random matrices, see the 
new review [7j reporting on this progress. For example, free probability 
has proven as an efficient tool for calculating the macroscopic level den¬ 
sity 0. With the help of orthogonal polynomials one could calculate alge¬ 
braic structures like determinants and Pfaffians, their kernels, and certain 
universal statistics on the local scale of the spectrum 0 m an E2 Eg. 
In particular products of random matrices drawn from Meijer G-ensembles 
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(the weight is essentially given by Meijer G-functions, see [14| for a defini¬ 
tion of these functions) exhibit a new kind of universal kernel in the hard 
scaling limit (microscopic limit around the origin). This limit is called 
Meijer G-kernel. Its name is reminiscent to the fact that the kernel essen¬ 
tially depends on Meijer G-functions. The “standard candles” of RMT, 
the Wishart-Laguerre ensemble [15] (xG/3E), the Cauchy-Lorentz ensem¬ 
ble [16] (£/3E), and the Jacobi ensemble [T7] (J/3E) are particular cases 
of Meijer G-ensembles. Also products of matrices drawn from these three 
ensembles are Meijer G-ensembles since this class of ensembles is expected 
to be closed under matrix products. 

Most results on the singular value statistics about product matrices 
are known for complex matrices (/3 = 2), only. The only exception, the 
macroscopic level density, can be computed for real (/3 = 1) and quaternion 
{(3 = 4) matrices with free probability [1] because they share the level 
density with (3 = 2 . However the local statistics of the singular values 
is still highly involved for (3 = 1,4 due to unknown group integrals like 
the Itzykson-Zuber integral [18j and its polynomial counterpart na mg. 
The projection formula recently proposed [20] circumvents such problems. 
This formula is a short cut of the supersymmetry method ED [22] and 
directly relates the original probability density with the weight in the dual 
superspace. 

After introducing the required notation in Sec. [2] we briefly review the 
projection formula in Sec. [3j Thereby we only consider the average of a 
product of characteristic polynomials to keep the calculation simple. We 
emphasize that the projection formula holds for all three Dyson indices 
/3 = 1,2,4 which is the strength of this approach. 

In Sec. 01 we demonstrate via the three ensembles, yG/3E, £/3E, and 
J/3E, how the projection formula works. Thereby we explicitly compute the 
well-known orthogonal polynomials for (3 = 2 and show that the average of 
one characteristic polynomial for (3 = 1 and the square root of a character¬ 
istic polynomial for (3 = 4 is apart from some shifts in the parameters the 
same as in the case (3 = 2 . Another example is presented in Sec. [5] where we 
generalize the approach to a product of independently distributed matrices. 
Also for product matrices we explicitly calculate the orthogonal polynomi¬ 
als in the case (3 = 2 . However the completely new results are the ones 
for (3 = 1,4 which are expressed in terms of integrals over Dyson’s circular 
ensembles (C/3E) [23] . In this way we show in Sec. [U] that the universality 
in the hard edge scaling limit holds for real and quaternion product ma¬ 
trices, too. We are also able to identify the non-linear c-models which are 
necessary when comparing the universal results with physical field theories. 


2. Preliminaries 


We consider rectangular random matrices which are either real (/3 = 1), 
complex {(3 = 2 ), or quaternion (/3 = 4). We are particularly interested in 
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the singular value statistics of a random matrix 

r M nx ( n+y ), j3 = i, 

W egl {p \n;n + iy) = ! c nx ( n+, '\ 0 = 2, (1) 

( e nx ( n+1/ ), 0 = 4 

distributed by P(WW^). We assume v = 0 in the following to keep 
the computations simple such that we choose the abbreviation gl^(n) = 
gl^(n;n). Nonetheless this restriction is not that strong since a product 
of rectangular matrices can be always rephrased to a product of square ma¬ 
trices p[2] . Examples of such induced measures resulting from rectangular 
matrices are given in Sec. [5j 

Since we choose the complex representation of the quaternion numbers 
H in terms of Pauli matrices we introduce the convenient parameters 


0 


4 

0 ’ 


7 = 


1, 0 = 1,2, ~ f 2, 0 = 1, 

2, 0 = 4, 7 \ 1, 0 = 2,4. 


( 2 ) 


For the sake of readability we restrict ourselves to partition functions of 
the form 


Z{M) = J d[W]P(WW^)deb 1 ^(WW i ® 1 ^ - M). 


(3) 


The fixed matrix M = {} has the dimension (777 x 771 ) < 8 ) (yk x 7 /c) = 
77 nk x 7777 .A;. It has to satisfy the symmetry 

m t = j In ® [t 2 ® lfe]Ml ra ® [r 2 <8) lfc], 0 = 1, / 4 x 

\ [r 2 ® In] ® lfe M [t 2 <8> In] ® lfc, 0 = 4, 

where r 2 is the second Pauli matrix. Other properties of M are not required. 

The partition function ([3]) needs an explanation. The determinant acts 
on the tensor space of (777 x 777 ) matrices containing the matrix WW^ and a 
space of dimension (ykxyk). In the case that M = l 7 n (g)diag (mi,..., m 7 fc) 
the determinant is a short hand notation for a product of characteristic 
polynomials of WW^ which is a well-known partition function in random 
matrix theory | 8 ]. The reason why we wrote this product in such an un¬ 
common, compact form is the application we aim at, namely the singular 

value statistics of matrix products. Then the matrix M does not take such 
a simple form. 

Another particularity of Eq. ({3]) which needs an explanation is the ex¬ 
ponent of the determinant, — 1 /( 77 ) an d the matrix dimensions. In the 
case of complex matrices (0 = 2 ), the exponent and the dimensions be¬ 
come self-explanatory since they become trivial, e.g. — 1 /( 77 ) = ~ 1 - 

When W is real (0 = 1) then WW^ is real symmetric and 77 x 77 dimen¬ 
sional. The space dual to the polynomials consists of self-dual matrices. 
The resulting Kramers degeneracy cancels the exponent 1/2 and doubles 
the dimension, k 2k. Exactly the opposite happens in the case of a 
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quaternion matrix W (0 = 4). Due to its quaternion structure the dimen¬ 
sion is doubled, n —>■ 2re. However the dual space consists of symmetric 
matrices. Since symmetric matrices may have also odd dimensions we do 
not need a doubling of the dimension k. The corresponding square roots of 
the characteristic polynomials are exact and, thus, a polynomial because 
the spectrum of WW^ is Kramers degenerate. Such a square root is known 
as quaternion determinant and is equivalent to a Pfafhan determinant [24j . 

An important ingredient needed for the supersymmetry method is the 
invariance of the probability density P under the transformation P{WW') = 
P(UWWW) for all U € U^(n) where 


U (/ 3 ) (n) 


O(n), 0 = 1, 

U(n), 0 = 2, 

USp(2n), 0 = 4. 


(5) 


Only due to this invariance it is possible to find an integral over a super¬ 
matrix whose dimension is independent of the ordinary dimension n and 
which yields exactly the same partition function as Eq. <[3j). This can be 
achieved in four steps which we briefly sketch in section [31 

For this purpose we have to introduce two supermatrix spaces and one 
ordinary matrix space. Let p,q,N S N, and U(p|g) and UOSp(p|2g) be 
the unitary and the unitary ortho-symplectic supergroup, respectively, see 
ESI [261 [27] . The space of rectangular supermatrices is defined by 


gl w (p|?;pV) = u (/3) (p + p r |g + g')/[ u(/3) (Pk) x u (/ 3 ) (pV)]> 


where u^(p|g) is the Lie superalgebra of the supergroup 


( 6 ) 


( UOSp(+)(p|2g), 0 = 1, 

U (/ 3 ) (p|g)={ U(p|g), 0 = 2 , (7) 

( UOSp ( - } (2p|g), 0 = 4. 


The coset is taken via the addition as a group action on the Lie superalge¬ 
bra. Therefore a matrix p e g\^(p\q;N) is ('ypl'yq) x ('yp\jq) dimensional 
and has the following form 


= Pbb P BF 
Pfb Pff 

The 7 px'yp dimensional boson-boson block pbb and the ygxyg dimensional 
fermion-fermion block ppp comprise commuting variables while the other 
two block contain anti-commuting ones. 

We employ the same notation for the two inequivalent fundamental 
representations of the supergroup UOSp(p|2g) as in [251 12ZJ where the su¬ 
perscripts indicate the transformation property under the complex conju¬ 
gation, i.e. 


P 


diag (l p , -it 2 <S> 1<j) pdiag (V,«t 2 ® l q >), 0 = 1, 
diag (~it 2 <8> 1 P , l q ) pdiag ( it 2 ® V, l q >), 0 = 4 


(9) 
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for p E gl^ 3 \p\q',p'\q') and 

rj* j diag (l p , —iT 2 <8> l g ) U diag (l p , zt 2 <8> 1 9 ), P = 1, nf) x 

\ diag (—zr 2 ® l p , 1 9 ) U diag (*t 2 <8> l p , l g ), P = 4 

for £/ E U^(p|g) C U( 7 p| 7 g). The two relations (0 and (fTUD are general¬ 
ization of the definitions of real and quaternion matrices to superspace. 
The ordinary matrix space announced is the coset 

f U(fc)/0(fc), p = 1, 

C/9E( 7 fc) = { [U(fc) x U(fc)]/U(fc) ~ U(fc), P = 2, (11) 

U(2Jfe)/USp(2fc), £ = 4 


equipped with a normalized Haar measure dp(U ) induced by the Haar mea¬ 
sures on the defining groups. These three sets are the circular ensembles 
first studied by Dyson [23]. These co-sets are also the fermionic part of 
the supermatrices involved in the superbosonization formula [22]. Since 
we only discuss the average of products of determinants and not ratios 
superbosonization reduces to bosonization only involving the circular en¬ 
sembles m ■ Let us recall the properties of a matrix U E C/ 3 E( 7 &). The 
matrix U is unitary and satisfies the symmetries U T = U for P = 1 and 
U T = (r 2 <S> 1 fc)L r (r 2 <S> lfc) for P = 4. 

Also the superdeterminant and the supertrace play an important role 
in the ensuing calculations. They are defined via the ordinary determinant 
and trace and explicitly read 

det (pBB - PbfPffPfb) ^ ^ /irA 

Sdet p =---, Str p = tr pbb — tr/O ff (12) 

det pff 

for an arbitrary square supermatrix p E gl^ 3 \p\q', p\q) whose fernrion-fermion 
block pff is invertible. The definitions are chosen in such a way that 
many properties of the trace and the determinant carry over to super¬ 
space. For example the circularity Str AB = Str BA, the factorization 
Sdet AS = Sdet A Sdet B, and the relation In Sdet A = Str In A still hold 
for two arbitrary invertible square supermatrices A and B. The circularity 
property of the supertrace works for rectangular supermatrices, as well. A 
more profound introduction in supersymmetric analysis and algebra can be 
found in [28] . 


3. What is the Projection Formula? 

The projection formula in its general form projects functions living on 
a very large superspace to functions on a much smaller superspace [ 20 ] . 
In this way it directly relates the original weight P to a weight Q in the 
smaller superspace. Hence the projection formula is a short cut of the su¬ 
persymmetry method [20]. For our particular purposes the large superspace 
is gl^(n + 7 ZI 7 Z; n| 0 ) with l being an integer larger than or equal to kj 7 . 
The enlargement of the dimensions k 21 in the case k odd and P = 4 
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is crucial. The reason is a Cauchy-like integration theorem [29j E 3 first 
derived in a general framework by Wegner m which only applies to an 
even dimensional reduction of a matrix space in the case of (3 = 1,4. 

In the first step of deriving the projection formula we need the following 
version of this Cauchy-like theorem Enj 

p{ww f ) = — 73 ) 

f d[n]exp[-strnnt] 

with W € gl^(ra) = gl^ (n|0; n|0) and fl € gl^( 7 ZI 7 Z; n|0). The matrices 
are embedded as follows 


' w ' 


■ w' ' 

n 


n 1 


€ gl^(n + 7ZI7Z; n). 


(14) 


The second splitting in W' 6 gl^(n+ 7 /| 7 Z —fe; n|0) and fZ' e gl^(0|fc; ra|0) 
becomes relevant in the third step of the derivation of the projection for¬ 
mula. The measure d[!Z] is the product of all differentials of independent 
matrix entries of fZ. The normalization with a Gaussian is true because the 
proportionality constant is independent of P and thus can be fixed by any 
weight. 

In Eq. (11311 we have chosen a supersymmetric extension of P to the 
superspace gl^(n + 7 ZI 7 Z; n|0) which is by far not unique. However the 
final result is independent of this choice as already discussed in [12]. Such 
an extension indeed exists for a smooth distribution P. Since P is invari¬ 
ant under the group U^(n) we can apply the Cayley-Hamilton theorem 
implying that P can be expressed in matrix invariants like traces and de¬ 
terminants of WW^. Those invariants have invariant extensions, namely 
the supertrace and the superdeterminant, cf. Eq. (USD- 

ill the next step we rewrite the determinant in Eq. ([3]) as a Gaussian 
integral over a matrix V = {1 / a j} € gl(^(0|fc;n|0) which only consists of 
Grassmann (anti-commuting) variables [28] , 

det 1 / (^)(WIT t <S> lyfc - M ) (15) 

= / m exp[tr rirn/Tt _ M aVJ K.V t ‘] 

f d[V] exp[tr VV^] 

Then the partition function is up to a constant 


Z(M) oc 

with 

V = 


exp 


-yk 7 n 

-StiQQ^V- E E M a b,ijV a iV£j 

a, 6=1 i,i=l 

(16) 


0 0 
V 0 


v* = 


0 

0 0 


€ gl^)(n + 7 ZI 7 Z; n + 7 ZI 7 Z). (17) 
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The first ('yn + 277 1 — 7 k) rows and the last 277 1 columns of V are equal 
to 0. The change of the sign in front of the first term in the exponential 
function relates to the fact that Grassmann variables are anti-commuting. 

The integrals over V and fi can be interchanged such that we find the 
function 


P(V f V) = j d[Q]P(^ t )exp[-StrS20 t t> t y]. 


(18) 


The invariance of P(£l fit) = P(UQftfU^) for all U € (71 + 7^7/) carries 

over to a symmetry for P(V^V) = P{UV^VU^) for all U € U^^ro + yilyi). 
Therefore the following duality holds 


P(V*V) = P{VV f ) 


(19) 


which is the third important step of the derivation. Employing the defini¬ 
tion (1181) backwards the partition function is 


Z(M) oc J d[n]d[V}P(n&)ei 

oc j d[fi / ]Q(n , f2 ,t )det 1/ ^(l 7 n ®n / n ,t -M). 


7 k 7 n 

trWn'fi'V- Y Y M ot>P v a,v h ) 

a, b=l i,j =1 


( 20 ) 


In the last step we integrated over the remaining degrees of freedom W\ 
cf. the splitting m , which do not show up in the determinant. This 
integration yields the function 


oc 


d[W']P 


w'W t w'n't 

n'w'* fi'fi't 


( 21 ) 


This equation is the essence of the projection formula. The remaining 
things to do is cosmetics. 

We want to express the dyadic matrix fTfl't as a single square matrix 
U which is an element in C/ 3 E( 7 fc). Note that the circular ensemble really 
relates to the Dyson index /3 = 4//3 and not /? which originates from the 
symmetries fulfilled by V. 

Exactly this is done in the last step. We apply the superbosonization 
formula |22] which reduces to pure bosonization in our case. This yields 
the partition function 


Z{M) 


J d/i([/)Q(D)det 1 /^(l 7n <8> U - M)det~ n ^U. 


( 22 ) 


with the normalized distribution 


Q(u) 


fd[Wi]d[W 2 ]P 


W{W\ + w 2 w\ 
U l / 2 w\ 


W2U 1 / 2 

u 


J dfi(U)d[Wi]d[W 2 ] det”™/ 7 U exp[—Str (WiWf + W 2 W\) + tr U] 

(23) 
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The reduction of the integral (1211) to the final expression (1231) as an integral 
over the two matrices W\ € gh^(rt + 7 ZI 7 Z — k;n + 7 ZI 7 Z — k) and W 2 € 
gl^(n + 7 i| 7 Z — k: 0 |fc) was done in [27j and is skipped here due to the lack 
of space. 

We remark that apart from the case k odd and j3 = 4 the auxiliary 
parameter l can be chosen l = k/ 7 . Then the matrix W\ is an ordinary 
square matrix and W 2 is a rectangular matrix only consisting of Grassmann 
variables. 

4. Application to Standard Random Matrix Ensembles 

Three particular cases of Meijer G-ensembles are the Gaussian xG/3E, 
the heavy-tailed C/ 3E, and the compactly supported J/3E. We discuss them 
in subsections HU H21 and 14.31 respectively. These ensembles play im¬ 
portant roles in a vast of applications and cover a broad range of sys¬ 
tems usmainiisiiniiEi. 


4-1 . Wishart-Laguerre (Gaussian) Ensemble 
The first ensemble we consider is the xG/3E, 

Pwl(WW^) oc det U ^WW ] exp[-tr WW f /T 2 ] (24) 


with v £ No and T > 0. It is the oldest random matrix ensemble first stud¬ 
ied by Wishart m- The determinant in front of the Gaussian originates 
from a transformation of a rectangular matrix W' € gl^(n, n + u) to the 
square matrix W € gl^(n). Therefore one can understand Eq. (|24D as an 
induced measure H2J. The corresponding weight Qwl is given by Eq. ([23]) . 


Qwl{U) oc J d[Wi]d[W 2 ]Sdet^ 


W{W\ + W 2 W\ W 2 U G 2 


u x i 2 w\ 


u 


x exp[—Str {W{W\ + W 2 W\) + tr U/T 2 } 
oc det- v ^Ue tlU / r2 . 

Therefore the partition function ([3|) for Pwl(IEVE^) reads 


(25) 


ry .... fdii(U)det-( n+ "Vmdet l ^Xl in ®U-M)e trU / r2 

ZWL(M) = : - JmU) -' <6) 


This result agrees with the one derived in m- The normalization can be 
fixed by considering the expansion of the partition function for large M. 

The result (1261) exhibits nice implications. For example the case k = 7 
and M = ml^,, is equal to the orthogonal polynomials for (3 = 2 and 
to the skew-orthogonal polynomials of even order for (3 = 1,4, see [24] . 
Hence the contour for j3 = 2 is a representation of the modified Laguerre 
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polynomials L 


n , 


see [T^, i-e. 


y{5— 2 ,^— 1 ) 

Z WL 


(min) 


oc 


oc 



m) n e z/r2 



1 



These polynomials also appear for f3 = 1,4 if we set k = 1. Only the 
argument m is modified to 7 m. Interestingly, the case (3 = 4 is an average 
over a square root of a determinant which is equivalent to a Pfaffian. 

For the case k = 2y and M = l 7 2 7n <g) diag (mi, m 2 ) we find one of 
the kernels corresponding to the yG/3E [24]. When computing the contour 
integral (1261) we immediately find the corresponding Christoffel-Darboux 
formulas. 


4-2. Cauchy-Lorentz Ensemble 

The £/3E is the next case we want to study. It is defined by the proba¬ 
bility density [16], [ 20 ] 

Pcl(WW^) oc det J// ^VEIE t det _/i (T 2 l 7n + WW f ) (28) 

with T > 1 , v G No and [i > k/ 7 + ( 2 n + i ')/7 — ( 77 — l )/2 for guaranteeing 
the convergence of the integral 014. It is a heavy-tailed distribution and 
was employed for modelling financial correlations [16]. 

The choice T > 1 is convenient for the projection formula but is not a re¬ 
striction at all because it only rescales the ensemble. The term det v ^WW^ 
can be again understood as a remnant of a rectangular matrix W' € 
gl^(n,n + v). However we underline that such a transformation from 
W' to W also changes the exponent fi. 

The weight for the dual space is calculated by Eq. (1231) . 


QcUU ) oc J d[Wi]d[W 2 ]Sdet^ 


W{W\ + w 2 w\ W2U 1 / 2 

U^wl u 


xSdet -'- 1 


r 2 l 7 n +7 7/|77^7fc + W M + W 2 W 2 ^U 1 ' 2 

u^wl r 2 i 7fc + u 

oc det _1/ ^[/det /i (r 2 l 7 j t + U) j c/[IEi]d[IE 2 ]Sdet V ^<W\W\ 

xSdet-** [r 2 i 7n+77W _ 7fc + WiW} + T 2 W 2 (r 2 l 7fc + U)~ l wl 
oc det- u ^Ude^- k ^- n ^(r 2 l^ k + U) (29) 


1 Note that the inequality satisfied by p in m contains a mistake which we have 

corrected here. The inequality can be found by performing a singular value decom¬ 
position of W and then reading off the algebraic behaviour at infinity. 
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In the last step we have rescaled W 2 —> W 2 (r 2 l 7 fc + U) 1 / 2 such that the 
remaining integrals are independent of U. Thereby we recall that the 
Berezinian (Jacobian in superspace) is (V 2 l^k + U). Hence 

we end up with the partition function 

„ n-\-u k n _ 1 

^ f dfj,(U) det t Udet^ t 7 (T 2 l^ fc + t/)detrr (l 7n <g) U — M) 

CL ~~ /d^Od^P^^ ' 

(30) 

Starting from this formula one can again easily deduce the orthogonal or 
skew-orthogonal polynomials, the kernel involving two characteristic poly¬ 
nomials, and the Christoffel-Darboux formula associated to this kernel. 
For example the orthogonal polynomials corresponding to the complex 
C/3E = LUE is 

oc j> dzz~( n+u+1 \r 2 + zy- n - 1 {z-m) n (31) 

n 

E l / m\J 

J=u jK n ~ j)\(y + j)\T[fi - n - u - j] \ T 2 ) 

This polynomial can be understood as a Jacobi polynomial when analyti¬ 
cally continuing the parameters to negative values, cf. Eq. m- The same 
polynomials pop up for /? = 1,4 when setting k = 1. This time we have 
only to change the exponent // —>• 7// — 7/7 + I. 


4-3. Jacobi (Truncated Unitary) Ensemble 

The J/3E is defined by [l7j 

Pj(IFVF t ) oc det^kFJF t det K (T 2 l in - W , W t )0(T 2 l 7n - W f ), (32) 

where u G No, ac > — 1 /( 27 ). The Heaviside step function for matrices 0 
is unity if the matrix is positive definite and otherwise vanishes. Again 
the scaling T > 1 is only introduced to avoid problems with the contour 
integrals in the dual space. In the case 77 /r G No the random matrix W 
distributed by Eq. (1321) can be understood as a truncation of an orthogonal 
(/5 = 1), a unitary (/3 = 2 ), or a unitary symplectic (j3 = 4) matrix, 
respectively, see mm 

To apply the projection formula we have first to find the supersymmetric 
generalization of the Heaviside step function. For this reason we write this 
function as 0(T 2 1 7 „ — W^W). Then it is clear that this function reads in 
terms of the supermatrix H as ©(r^l-„ — because the dyadic matrix 
has still an ordinary dimension and can be embedded in the space of 
7 n x 7 n matrices by a Taylor expansion in the Grassmann valued matrix 
entries. Such a Taylor expansion is always finite since Grassmann variables 
are nilpotent. Hence we do not have to fear any problems of convergence. 
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Let n,p,q € N and V € gl^(p|g; n|0). Then the extension of the 
Heaviside step function is done by a limit, 


0(l 7n - TV) 



(33) 


This limit vanishes if one or more eigenvalues of the numerical part of 
the dyadic matrix V^V is larger than 1. We emphasize that indeed only 
the numerical part matters and not the nilpotent terms because of the 
Taylor expansion in the latter. In the next step we employ the duality 
tr e^ tV<< ' = 7(77. — p) + yg + Str e^ eVyt . We have 

0(U-LV) = lim(l + e- e )^- n )-^Sdet- 1 (l 7P | 7 „ + e- e e eyyt ) 

= lim det _1 (l 7p i 79 + e~ e {e tVy] }bb) 

= 0(l 7p -{HW } nu m) _ (34) 

The Heaviside step function is only taken for the numerical part 
of the boson-boson block of the dyadic matrix VV^. Any expansion in 
the nilpotent terms yields a polynomial in e which are suppressed by the 
exponential e~ e . This implies that the other three blocks of the supermatrix 
e eVy1f canno t contribute because they are polynomials in e. The boson- 
boson block is {e tVyt }BB = ' (1 + /(e)) with / a polynomial with 

/(0) = 0. Therefore Eq. (1341) is the correct generalization of the Heaviside 
step function to the superspace. Interestingly the Taylor expansion in the 
nilpotent terms have no influence on the Heaviside step function. But 
this behaviour has to be expected because the Taylor expansion can only 
have an effect on the boundary. Only there one or more eigenvalues of the 
numerical part are equal to 1 where the value of the function may 

change. However, the supersymmetric Heaviside step function vanishes at 
the boundary, too, due to the expansion in the nilpotent terms yielding an 

inverted polynomial in e, e.g. det -1 (l 7J ,i 7l? + e~ e {e eVV ^}BB) ^ 1P 

1 //(e) 0 with / a polynomial. 

We employ Eq. ()34[) in our setting and recognize that the matrix U 
is not a part of the boson-boson block of the matrix argument of Pj in 
Eq. (1231) . Hence, the function in the dual space is 


Qj(U) oc J d[W 1 ]d[W 2 ]Sdet v ^ 


W X W\ + W 2 W ] 2 W 2 U 1 / 2 

U x ! 2 w\ u 


x Sdet K 


r 2 i 7n+77 * l77 *_ 7fc - w^wl - w 2 wl 

U l / 2 wl 


X0(r 2 l 7n+77Z - {IT, w/ + w 2 wl }| U B m ) 
oc det- l// ^Hdet- K - fc/7 - n/ ^(r 2 l 7fc - U ). 


W 2 U 112 

r 2 i 7 fc - u 


( 35 ) 
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We underline that the boson-boson block of W 2 Ulf only consists of nilpotent 
parts such that it does not contribute to the Heaviside step function. The 
corresponding partition function is 


n-\-u k n _ 1 

f d/i(U) det t Pdet 1 ^ {T 2 l^ k ~ U)det^ {l in ®U — M) 

J _ Jdfi(U)det~^ n+ ^^Udet~ K - k ^- n ^{T 2 l^ k - U) ' 

(36) 

One can readily check the correctness of this result by calculating the or¬ 
thogonal or skew-orthogonal polynomials and the kernel involving two char¬ 
acteristic polynomials. For example, with the help of the residue theorem 
we generate the polynomials, 


<* (37) 


oc 


y- ^[ri + K + u + j + 1 ] /my' ((tiI/) /2m _ \ 

e n jKn-j)\(u + j)\ \ r 2 ) n Vr 2 J ’ 


where P^' u ^ are the Jacobi polynomials with respect to the weight (1 — 
x) K (l + x)"@(l — x 2 ), see [13]. As in the case of the £/3E we find the 
same polynomials for f3 = 1,4 and k = 1 when replacing the exponent 
K —>■ 7 K + 7/7 — 1. 

We also obtain the well-known Christoffel-Darboux formula of the Ja¬ 
cobi polynomials by setting k = 2, /3 = 2, and M = l n <g> diag (mi, m 2 ). 
Then the integral reduces to a double contour integral after diagonalizing 
U. 


5. Application to Product Matrices 

The computation of the partition function for a product of L matri¬ 
ces W —> = nt=i Wj = Wi''' kHz, independently distributed by 

P(ITVF^) —> riy=i Pj{WjWj) works in a similar way as for a single matrix. 
Starting from the partition function 


Z n (M) 



det 1/(77) 


w {l) (w^y ®i^ k 



det k hw^ L ~ v > (w^- 1 ^ (38) 


xdet 1/(77) 


W L W[ ® U k - X-^MY-\ 


with Xl~i = W^ L ^ (g) l^ k and Yl~ 1 = (VF^ L ^)t (8) lyfc, we apply the 
projection formula for Wl after replacing the matrix M —>• XJ^_ j MY^ { . 
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Then we obtain 


Z U (M) 



d^(U L )Q L (U L )det~ n ^U L 


xdet 1/(77) 




(8 )U L - M 



d^(U L )Q L (U L )det k ^W^ 



t 


1 

x det 


Wl-xW]^ ® l, fk - Xl\MY~\ 


(39) 


where Xl ~2 = W^ L ~ 2 ^ <g) \fTJp, Yl -2 = (fF( i_2 ))t ® \fUp , and Ql is com¬ 
puted as in the projection formula (1231) . This procedure yields a recursion 
resulting in the following expression for the partition function, 


Zu(M) 



det ~ n ^U L • ■ ■ U\ 


xdet 1 ^ 


l 7 n ® VTTl • • • y/lhUiy/lh- ■ ■ - M 


m 


where each matrix Uj is an element in the circular ensemble C/3E(7 k). 

In the final step we replace C/j = \jU k ■ ■ ■ yjUj + \Uj y/Uj + \ ■ ■ ■ \fTJp 
which preserves the symmetries such that t/j € C/3E(7 k). For this pur¬ 
pose we use two facts. First, the Haar measure is invariant under dfr(U) = 
dn(VUV T ) for all V € U( 7 k) resulting from the fact that the explicit 
form of the Haar measure of CfdE^k) is d/j,(U) oc de\F k / 1 ~^~~d/ 2 ud\U] 
with d[U] the product of the differentials of all independent matrix en¬ 
tries [23 22]. Second, the weights Qj are also invariant under Qj(U) = 
Qj(VUV^) for all V £ U ^\k). Hence these weights have an expression 
in terms of functions of matrix invariants. With the help of a slight abuse 
of notation one can say that the weights Qj satisfy a cyclic permutation 
symmetry, Qj(AB) = Qj(BA) for any two matrices A,B € U( 7 &;). 

Finally, we find the result 


Z U (M) 



det" n/7 H(det 1/(77) [l 7n <g> U[ — M] 

(41) 


with U' L+1 = l 7 fc. This result is surprisingly compact. It also reflects the 
nature of the original product of matrices which is equivalent to a Mellin-like 
convolution in a matrix space. Also the dual space exhibits this structure 
of a Mellin-like convolution. 

As an example we calculate the orthogonal polynomials {k = 1) of a 
product of l complex xG/3E = xGUE, Eq. (f24j) . Lcl complex Cj3E = 
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LUE, Eq. d2U), and Lj complex J/3E = JUE, (124|) . We assume this product 
to be ordered, i.e. first the Wishart-Laguerre, then the Cauchy-Lorentz, and 
finally the Jacobi matrices. The result does not depend on this ordering, 
see the discussion in [12| . Then the orthogonal polynomials are 


z w.cLij( mi ")«i n t 2 


dz* 


-1 


. J = 1 


z 3 +1 


n v,- 1 

r? *j+i 

e j 


7 ? J 


(42) 


/iwL+^CL 


dZj_ 

"u+i" 

v i 

n 

j. to 

+ 

d? 

1 

Z 3 

L Z 3 \ 


L *7+1J 


1 


oc 


cfe 9 - 


x 




( ■^WL+^CL+^J 

n . 

j'=Lwl+^CL+1 J 

n B r[n + «a + i' + j + 1 ] 


' z 3+ 1" 

V 3 

1 

1 

CN v 

t-H 

l Z 3 \ 


L u+iJ 


- K '-r 


my 


J -^i ! ( n -j) ! (n„(^a+ j) ! )(ria r [Ma-n-I/- j]) V T 2 / 


with z Lwl+Lcl+Lj+ i = 1, n'j = Hj - n, «'■ = Kj + n, and T 2 = [], T 2 . The 
product of the Gamma functions runs over the possible values for u a , K a , 
and p a . The polynomial (1421) is a hyper geometric function and, thus, a 
Meijer G-function m It agrees for certain values of the parameters TwLj 
Tcl> an d Lj with known results mum What is completely new are the 
results for /3 = 1,4 and k = 1 which are essentially the same polynomials. 
Here the other approaches failed because of unknown group integrals. 


6. Hard Edge Scaling limit of Product Matrices 

Up to now every calculation was done for finite p such that we made 
no approximation and the projection formula was exact. However to make 
contact to physical systems and universality we have to zoom onto the local 
scale somewhere of the spectrum. A very prominent scaling is the one to a 
vicinity around the origin also known as the hard edge scaling limit. 

As a simple but non-trivial example, we choose the matrix product of 
the previous section with the source M = <g> rh/\n(l\a(Va - 

n h))^Wai K a + n/ 7 ))]. In particular we consider the scaling limit n —>• 00 
and Vj, Jij = {pj/n — 1 / 7 ), Kj = ( Kj/n + 1 / 7 ), and m fixed. Then one can 
easily show that the asymptotics of each weight, regardless what kind of 
random matrix we consider, is 

Qj ( aU) "(x 1 det^Ue* u (43) 

with a = T 2 for xG/3E, a = T'j/(npj) for £/3E, and a = T 2 -/(rikj) for £/3E. 
After a proper rescaling of the matrices Uj the partition function (1411) takes 
the asymptotic form 

J\dpL(U'j)det^^Uj \ e trUL +^=i tTU i u i+i~ tr ac7 i _1 (44) 


Zn(M) n « 1 J 
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with L = Lwl+Ecl+Ej. We underline that no saddlepoint approximation 
is needed for this limit. Hence the matrices Uj are still elements of the 
circular ensemble C/3E(7fc). 

For /3 = 2 the partition function (1441) yields the Meijer G-kernel of a 
product of matrices drawn from yG/3Es, cf. m- This can be seen by di¬ 
agonalizing the unitary matrices, applying the Itzykson-Zuber integral [,18j| 
and finally integrating over a determinantal point process. The entries of 
the resulting determinant are Meijer G-functions. Our result emphasizes 
the conjecture that also this kernel is universal. Indeed we could also have 
chosen another scaling which still leads to a hard edge scaling limit. Then 
we would get finite rank deformations of the result (I44[) which was recently 
discovered for a product of truncated unitary matrices in m- Nevertheless 
the limiting kernel is still a Meijer G-kernel but with other parameters. 

From a physical point of view can one ask for the non-linear u-model 
corresponding to the partition function (1441) . In this framework the func¬ 
tion in the exponential function is identified as the potential. The inte¬ 
gration domain C/3F, L (~/k) is the coset of the “flavour” group which keeps 
the “massless Lagrangian” (m = 0) in the full theory at finite “volume” 
n invariant divided by the group which keeps the ground state invariant. 
As in the case L = 1 the theory is spontaneously broken. For a product 
matrix the “flavour” symmetry at finite “volume” n is \] L (^/k) for /3 = 1,4 
and [U(fc) x U (k)\ L for (3 = 2 which can be readily checked by linearising 
the product in the matrices Wj. This group is spontaneously broken 

to \U^\k)] L and the source term for its condensate is the “mass” fh. This 
non-linear a model generalizes the one for the Wishart-Laguerre ensemble 
which were found in QCD m and mesoscopic systems [25! . Especially the 
coupling between different Uj is reminiscent but not completely the same 
as the recently proposed chiral Lagrangian for high density QCD [32] . 

7. Conclusions 

We briefly presented the projection formula [20] for averages over prod¬ 
ucts of characteristic polynomials which is a short cut of the supersymmetry 
method mm- The general results found by this approach were demon¬ 
strated in the case of Wishart-Laguerre (xG/3E), Cauchy-Lorentz (C/3 E), 
and Jacobi (J/3E) ensembles, in particular we rederived the corresponding 
orthogonal polynomials for (3 = 2. These polynomials are essentially the 
same when averaging over one characteristic polynomial for /3 = 1 and over 
a square root of a characteristic polynomial for /3 = 4. 

Moreover we generalized the projection formula to products of matrices. 
Since the projection formula works in a unifying way for all three Dyson 
indices (3 = 1,2,4 this approach is an ideal alternative compared to other 
methods like orthogonal polynomials and free probability when studying 
real or quaternion matrices. Note that up to now free probability only 
applies to global spectral properties and to use orthogonal polynomials we 
need to know group integrals like the Itzykson-Zuber integral m or its 
polynomial counterpart muni. The projection formula circumvents this 
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problem. In particular we were able to show that the spectral statistics 
at the hard edge are the same for products of completely different random 
matrices only depending on the number of matrices defined and their indices 
,ui encoding the rectangularity of the matrices. This was done for all 
three cases /3 = 1,2,4 and underlines the strength of the projection formula 
where other methods fail. In the complex case (/3 = 2) we easily deduce 
from our results those for the Meijer G-ensembles studied in mum 

The projection formula also enabled us to identify the non-linear a- 
models and the symmetry breaking pattern for product matrices and de¬ 
rived the potential of the Goldstone manifold. This result is completely 
new and shows what the effective theory associated to such a product ma¬ 
trix would look like. In particular one can understand a product matrix by 
itself as a discrete one-dimensional system. Therefore our results show one 
way to generalize the zero-dimensional RMT to a one-dimensional theory. 
Indeed via the DMPK equation such a link to a one-dimensional system 
was established [3j, though there a different limit is considered. 
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